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ABSTRACT 

The neurons work together to process the information in nervous system that in 
turn controls the animal behaviors. This is an attempt to understand the dynamics 
of one of the basic components of central nervous system i.e., neuron. Exploring 
the dynamical characteristics of neuron is very helpful in uncovering the mechanism 
behind brain functions. Along with this, revealing the parallel processing in brain 
requires knowledge of coupling and synchronization among different neurons. In this 
thesis work, the dynamical features and effect of coupling strengths on synchronization 
in network of neurons are studied. Bifurcation analysis are carried out to find out the 
variation in excitability considering stimulus intensity as the bifurcation parameter. 
Besides the dynamics of point neuron, the signal propagation in the neuron structure 
is also investigated. The neuron structure consists of amalgamation of dendrites, soma, 
axon and synapse each exhibiting different dynamics. The study of signal propagation 
in complete neuron structure is helpful to depict the dynamics of biophysical activities. 
The effect of random opening and closing of ion gates and ephaptic interaction are taken 
care of by adding white Gaussian noise in the analysis. A different method to evaluate 
signal propagation is proposed which is based on calculation of Fourier coefficients from 
the generated wave. 
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Notation 

Description 

Dimension 

t 

Time 

rns 

T 

Time Constant 

ms 

X 

Distance 

mm, cm 

v m 

Membrane Potential, Action Potential 

mV 

Vte* 

Resting Membrane Potential 

mV 

R, R,n 

Membrane Resistance 

Kil - an 2 

Cm 

Membrane Conductance 

mS/crn 2 

Rnj 

Stimulus 

nA 

k B 

Boltzmann constant 

m 2 Kg/(s 2 K) 

T 

Absolute Temperature 

K 

<1 

Charge of a single Proton 

C 

7^ 

O* channel 

Maximum conductance of an Ion Channel 

mS/crn 2 

9 channel 

Conductance of an Ion Channel 

mS/crn 2 

Z 

Valence of Ion 


[B] 

Concentration of B ions 

m.M 

Pa 

Activation Variable 


Pda 

Deactivation Variable 
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Notation 

Description 

Dimension 

G 

Number of Activation Channel 


C la 

Number of Deactivation Channel 


(*Na 

Maximum Conductance of Sodium Channel 

mSf cm 2 

C K 

Maximum Conductance of Potassium Channel 

mS/crri 2 

G Ca 

Maximum Conductance of Calcium Channel 

rnS/cm 2 

C L 

Leakage Conductance 

mS/crri 2 

On a 

Conductance of Sodium Channel 

rnS/cm 2 

Ok 

Conductance of Potassium Channel 

mS/cm 2 

(JCa 

Conductance of Calcium Channel 

rnS/ cm 2 

m 

Activation Variable for Sodium Channel 


h 

Deactivation Variable for Sodium Channel 


77 ,, W 

Activation Variable for Potassium Channel 


tllpo 

Steady state Calcium Activation 


Woo 

Steady state Potassium Activation 


Tw 

Time Constant for Morris-Lecar Model 



Rate Constant: Transition from closed to open state 

s~ l 

p. 

Rate Constant: Transition from open to closed state 

s- 1 

El 

Equilibrium Leakage Potential 

mV 

Enci 

Equilibrium Sodium Potential 

rnV 

E k 

Equilibrium Potassium Potenti 

mV 

Eca 

Equilibrium Calcium Potential 

mV 

E 

^ syn 

Synaptic Potential 

mV 

Qsyn 

Synaptic Conductance 

mS/cm 2 

t'peak 

Time to reach peak value for NMDA Current 

ms 



Notation 

Description 

Dimension 

Tl, T‘2 

Time Constants: Non-NMDA Synapse 

ms 

<jc' j 

Coupling Conductance between two Neurons 

rnS/cm 2 

r 

Equilibrium point value of any Variable i /> 

Same as 

j 

Jacobian 


A 

Eigenvalue 


1'L 

Intracellular Resistivity 

KD. — mm 

hn 

Membrane Current 

riA 

a 

Radius of Cable 

/mi 


Space Constant 

mm 

*k 

Longitudinal Resistance 

KQ 

fi 

Compartment N umber 


g^' 

Coupling Conductance between two Compartments 

mS/cm . 2 

V 

White Gaussian Noise 

dBW 

«fc, h 

Fourier Coefficients 


c 

Speed of Wave 

an /ms 

VD 

Potential of Dendrite 

mV 

I< 

Constant 

mS 

RD 

Coupling Resistance between Soma and Dendrite 

KQ 

Rs 

Soma Resistance 

KQ 

V s 

Summed Potential at Soma 

mV 
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Chapter 1 


Introduction 


The study of dynamical behavior of the nervous system has drawn huge attention 
during the last few decades. The main aim of this study is to find out the functioning 
of brain. Despite a lot of attempts, the proper mechanism of the brain is still unclear. 
The smallest integral part of any nervous system is neuron. Exploring the behavior 
of neuron will be very helpful to uncover the details of the brain dynamics. Brain 
dynamics include information processing in neurons such as generation of spike, its 
propagation, communication among neurons, bifurcation analysis, presence of chaos, 
synchronous behavior etc. In this thesis, effort is made to analyze the synchronous 
behavior of neuron and information transmission in nerve cells by applying nonlinear 
dynamical analysis techniques. 

1.1 Motivation 

The brain allows the synchronous functioning of neurons so that the neuronal system 
work more efficiently. Present trends of research include investigation of the behavior 
of neurons when considered in a network. Researchers have concentrated their study 
mainly on the dynamics of two coupled neurons. The extension of this study on 
increased number of neurons may yield significant results. Also, the studies have been 
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carried out on different parts of neurons viz. dendrites, soma and axon separately which 
helps in knowing the flow of information from one end to the other. The response for 
complete neuron structure is required in order to understand the phenomena occurring 
in small networks of nervous system. A better approach to model the signal propagation 
is needed that would be easier to compute when compared to previous methods. 

1.2 Literature Survey 

The simplest model of the nervous system is given by the dynamics of a point neu- 
ron. The first attempt, to predict the neuronal behavior was given by Lapicque in 
1907 [22]. He modeled the membrane as an RC circuit and introduced a. spike arti- 
ficially if the membrane potential crosses a threshold value. The model is known as 
Integrate- and-Fire model. The effect of ion channels are not considered in this model. 
In 1952, Hodgkin and Huxley proposed a biophysical model of the active membrane 
that incorporates the dynamics of ion channels [17]. Hodgkin- Huxley model considers 
the membrane current as a composition of leakage current, sodium current and potas- 
sium current. This model is four dimensional and hence complex and computationally 
inefficient. FitzHugh [12] proposed a simplified version of the Hodgkin-Huxley model. 
Nagumo et. al. derived the equations those qualitatively describe the events occurring 
in an excitable neuron [27]. The models proposed by them are combined together to 
give a two dimensional model known as FitzHigh-Nagrmo model. This model repli- 
cates the biophysical activities in terms of fast and slow variables. FitzHugh-Nagumo 
neuronal model is very simple and widely used by different researchers to understand 
the brain phenomena. Because this model does not represent any physical membrane, 
Morris and Lecar proposed a two dimensional model that portraits the dynamics of 
barnacle muscle fiber [25]. This model, known as Moms-Lecar model, considers the 
effect of calcium channel. It has many variants depending upon its parameter value, 
the values given by Rinzel [35] are widely used for analysis. All these models represent 
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either Iype-I or lype-II membrane classified based on their firing frequency. A vari- 
ant of neurons also exhibits the bursting phenomenon that is not captured by these 
popular neuron models. Rinzel proposed a modified three dimensional version [33] 
of FitzHugh-Ngumo model known as FitzHugh- Rinzel or Modified FitzHugh-Nagumo 
model. It shows the bursting in neurons. 

For parallel processing in brain, the neurons must synchronize with each other. 
The synchronization allows efficient execution of any task. There have been physi- 
cal evidences of coupling in nervous systems. A study of coupling in the developing 
mammalian retina has been carried out by Penn [28], Gerloff et. al. showed that 
task-related coherence reflects inter-regional functional coupling of oscillatory neuronal 
activity [14]. The recent trends in computational neuroscience research is to study the 
synchronizing behavior among neurons. An experimental study on a pair of biological 
neuron that interact through naturally occurring electrical coupling shows that the 
natural coupling synchronizes slow oscillations of membrane potentials [10], Ucar et. 
al. have studied the effect of coupling in FitzHugh-Nagumo system and designed an 
active controller that synchronizes the weakly coupled neurons [42], 

The local neural dynamics is determined by voltage-and ligand-gated ion channels 
and feedback between densely interconnected excitatory and inhibitory neurons. With 
modulation of long-range synaptic coupling, the system undergoes a transition from 
independent oscillations to stable chaotic synchronization [4]. Bifurcation analysis pre- 
dicts the effect of parameter variation in any dynamical systems [38]. Rinzel and Er- 
mcntrout studied the bifurcations in the Morris-Lecar model by treating the externally 
applied dc current as bifurcation parameter [34]. A bifurcation analysis and investiga,- 
tiou of presence of chaos is studied in [13], [19] and [24], Like single neuron models, 
the bifurcation analysis on coupled neuron models also reveals significant informations. 
In [2], a study is carried out in excitable chaos in coupled FitzHugh-Nagumo model. 
Kavasseri has investigated bifurcations in synaptically coupled Morris-Lecar model [21]. 
Similar to chaos and bifurcation, the study of noise is also important in nervous systems 
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as there are many sources of noise such as Synaptic bombardment , Ephaptic Interaction 
and uneven opening and closing of ion channels. Noise may change the behavior of the 
system drastically. The stochastic resonance helps in signal detection of neurons [37] . 
Tuckwell and Rodriguez showed both analytical and simulation results for stochastic 
FitzHugh-Nagumo model [41]. The authors in [32] carried out a similar analysis on 
Hindmarsh-Rose model. 

All of these studies have been carried out on point neuron models. However, there 
is a need to carry out analysis on axon and dendritic tree. The neuronal spikes show 
different characteristics when they travel along active axon or passive dendritic trees. 
There is a delay as well as decay in spike when it propagates through passive mem- 
brane [9]. The delay depends upon the speed of the spike which in turn depends on 
the topology of the cable. The speed of action potential in a giant axon of squid is 
18.8 m/s [26]. The axon stem brandies out near soma. This results in a change of 
speed. The diameters of the daughter branches do follow a pattern. A method to find 
out these diameters is given in [6]. The ion channels embedded in axon helps in proper 
propagation of action potential. The range of ion concentration that should be present 
for a successful propagation of signal is studied in [43]. Gordon and Welsh [15] dis- 
cussed about ion channels embedded in axon. The axons of different neurons dissipate 
a certain energy to produce a spike. Authors in [3] studied the energy requirement 
for spike propagation. The Ephaptic Interaction [18] depolarizes the axon and thus 
helps in signal propagation. This interaction occurs between two axons closely located 
and touching each other. The interaction between nerve fibers are also studied in [7] 
and [16]. Arvanitaki suggested that Ephaptic Interactions presume electrical fields in 
extracellular space generated by active neurons induce current flow in adjacent neu- 
rons [1]. There may be generation of spike in the dendrites if it is spine studded. An 
analysis supporting this phenomenon has been carried out by Coombes and Lord [23] . 
The traveling waves in dendrites may also be periodic [8]. The periodic response in 
axonal propagation has been found by Terakawa [40] 



1.3 Organization 


The thesis work has been organized in four chapters. 

Chapter 2 discusses on neuron and its dynamical behaviors. It presents a brief 
introduction on neuron structures and describes its characteristics. The dynamics of 
active and passive membranes and that of synapse has been studied. The effect of cou- 
pling on neuronal behavior is found out by investigating the synchronizing conditioxrs 
and bifurcation analysis. 

Chapter 3 covers the studies on the signal propagation in neuron. It focuses the 
important aspects of the models for signal propagation in dendrites and axon. Effect 
of noise on axon has been examined. A new method based on Fourier series has been 
proposed to model the axon. Finally, the picture of signal propagation and neuronal 
communication have been portrayed with the help of complete neuron structure. 

Chapter 4 gives the conclusions based on findings and suggests some future scopes. 



Chapter 2 


Neuron and Its Dynamics 

2.1 Structure of Neuron 

The information processing in nervous system is perhaps the most existing complex 
phenomenon in the world. The basic building block of the nervous system is neuron. 
According to biologists, a human brain consists of approximately 100 billion neurons 
with each neuron having at least 10,000 connections with other neurons. The communi- 
cation between two neurons usually takes place through synapse, a gap separating two 
neighboring neurons. Sometimes the neuron communicates through Ephaptic Trans- 
mission which is nonsynaptic electrical interaction [22]. The neurons are similar to the 
other cells in the bodies [46] as 

• Neurons are surrounded by a cell membrane 

• Neurons have a nucleus that contains genes. 

• Neurons contain cytoplasm, mitochondria and other organelles 

• Neurons carry out basic cellular processes such as protein synthesis and energy 
production. 


6 
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On the other hand, neurons have some properties that distinguish it from other cells 
and these properties are 

• Neurons have specialized extensions called dendrites and axons. Dendrites bring 
information to the cell body and axons take information away from the cell body. 

• Neurons communicate with each other through an electrochemical process. 

• Neurons contain some specialized structures eg. synapse and chemicals eg. neu- 
rotransmitters. 



Figure 2.1: Architecture of a neuron 


A neuron has three major morphologically defined portions: soma, dendrite and 
axon [20]. Each of these parts makes a specific contribution to the processing of nerve 
signals. Figure 2.1 1 shows the structure of a typical neuron found in vertebrates. The 
soma also known as cell body consists of the cell nucleus and perikaryon. It is usually 
50 /mi or more in diameter. The dendrites branch out from soma in treelike fashion. 
These are short nerve fibers and highly branched. The dendrites exist in many per cell 


‘Figure 2.1is taken from [47] 
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and covers approximately 80 — 90% of the total surface area of neuron. Axon is slender 
and may Ire very long. It has a diameter ranging from 0.2 to 20 pm and a length up to 
lor. It has several branches at its end near synapse. Sometimes, the axon is covered 
with myelin sheath consists of fat-containing cells that insulate the axon. The gaps 
(approximately 1 pm wide) formed between myelin sheath cells long the axons are 
called Nodes of Ranvier. A neuron has only one axon. It begins at the Axon hillock. 
The dendrites rocicve inputs from the axon terminal of other neurons at synapse. The 
electrical signals travel along the dendrites to soma and aggregation takes place there. 
The summed signal travels down the axon to synapse where it stimulates the dendrites 
of other neurons. The myelin sheath allows rapid and efficient transmission of impulses 
along the nerve cells. 

Cell morphology may be classified according to the number of elements emanat- 
ing from the cell body [36] as unipolar cells, bipolar cells and multipolar cells. The 
unipolar cells are generally found in invertebrates while multipolar cells are dominant 
in vertebrate nervous systems. Neurons can be further classified by the direction that 
they send information: 

• Sensory or afferent neurons that send information from sensory receptors viz. 
eyes, ears, nose, tongue and skin toward the central nervous system. 

• Motor or efferent neurons send information away from the central nervous system 
to muscle or glands. 

• Interneurons send information between sensory neurons and motor neurons. 

Each type of neuron poses different characteristics. The neurons are mainly classified 
as Type.-I , Type-II, and Type-Ill based on their dynamical response. Among these, the 
dynamics of Type-I and Type-II neuron models have been considered for the study as 
the most of the neurons follow either of the two characteristics. 
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2.2 Neuron Modeling 

2.2.1 Membrane Equation 

The cell membrane is a lipid bilayer, usually 3-4 ran thick, that is essentially imper- 
meable to most charged molecules. A cubic micron of cytoplasm might contain 10 10 
water molecules, 10 s ions, 10 7 small molecules and 10 5 proteins [9]. The mobility 
and repelling nature of the ions and the molecules that carry charges may result in 
a build up of negative charges inside the membrane. Electrostatic forces attract an 
equal density of positive ions from the extracellular- medium to the outside surface of 
the membrane. The insulating nature causes the membrane to act like a capacitor 
by separating the charges lying along its interior and exterior surface. Numerous ion 
channels are embedded in the cell membrane. The ion channels are highly selective 
in most of the cases allowing only a single type of ion to pass through them. These 
ion channels eventually lower the effective membrane resistance to a value about 10000 
times smaller than that of a pure lipid bilayer. When the neuron is inactive, the excess 
internal negative charges cause the potential inside the cell membrane to be negative. 
This potential is known as Resting Potential and denoted by Vr (M . 



V 

Rest 


TrFn 

Figure 2.2: Circuit representation of a passive membrane 


Thus a membrane can be represented in most simple way by a parallel combination 
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of a capacitance with a series connection between a resistor and a battery as shown in 
Figure 2.2. The dynamics of this linear model can be found out by applying KCL at 
node A. It results in Equation 2.1 

= + V Resl + ( 2 - 1 ) 

Where r = RC. The solution of Equation 2.1 for a step current input with magnitude 
/o is given by Equation 2.2 

V m (t) = RI 0 (l - e-^) + Vjtut (2-2) 



Figure 2.3: Response of passive membrane 

If the current injection is stopped at time t. 0 ff after the membrane potential reaches the 
steady state value, then the membrane potential returns to V Rf , at with an exponential 
time course given by, 

V m (t) = Rio + V Rest (2.3) 

The response of the passive membrane modeled using a simple R-C circuit is shown in 
Figure 2.3. The red, blue and black lines correspond to membrane potential at different 
time constants n, n, and r 3 respectively with r 3 > r 2 > t x . The magnitude of injected 
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current is same for all the cases. It has been verified that, the membrane potential 
reaches its steady state value more rapidly for lower values of time constant. For lower 
values of time constants, the passive membrane of dendrites reacts to a synaptic current 
change more rapidly. The green line shows the response of the membrane potential 
with the time constant value T\ and subjected under inhibitory current. As depicted 
from Figure 2.3, the dynamics of passive membrane represents the subthreshold region. 

2.2.2 Point Neuron Models 

The first attempt to represent the neuronal activities in terms of mathematical equa- 
tions was formulated in 1907 by Lapicque. Different mathematical models have been 
developed later on to represent the biological activities in a more detailed manner. 
Out of these models, the models that are found more frequently in literature are 
Leaky Integrate and Fire (LIF) [22], Hodgkin-Huxley (HH) [17], FitzHugh-Nagumo 
(FHN) [12], [27] and Morris- Lecar (ML) [25] neuron model. The HH, FHN and ML 
models consider the dynamics of different ion gates. These ion gates are highly se- 
lective and usually allow only a single type of ion to pass through it. Each of these 
channels may be represented as a battery in series with nonlinear voltage dependent 
conductance as shown in Figure 2.4. The potential offered by the membrane may be 
found out through Nemst Equation given in Equation 2.4. 

_ Vt ( [ion(outside)} 
z 11 \ [i.on(inside)\ 

where z is the valence of ion and V T is typical membrane potential. It has a value 
between 24 mV to 27 mV for the typical temperatures of cold and warm blooded 
animals [9]. The membrane potential varies from -3 to +2 times of V r . 
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'flic expression for typical membrane potential is given by Equation 2.5 

V T = 2.5) 

q 

where ku is the Boltzrnan constant,, T is absolute temperature, and q is charge of a 
single proton. 

The channel conductance are given by the following expressions. 

(/channel ( I'm ) = Lr channel (p^Pda^ (2.6) 

where the channel has Q, activation and Q a deactivation ion gates. p a and p ( i a represents 
the probability that the corresponding type of ion gate is open at any time instant. 
These ion gates are voltage dependent and hence vary with any change in membrane 
potential. G c h<mnd is the maximum possible conductance of the channel. 

Among the above mentioned models, the HH neuron model represents the dynamics 
of Type-II and ML model represents the same of Type-I neuron model. Therefore these 
two models are considered for analysis to cover up the different types of behavior 
observed in a neuron. 

2.2.3 Hodgkin-Huxley Model 

Hodgkin and Huxley proposed a neuron model [17] from the experimental data obtained 
by artificially stimulating the giant axon of squid. This model considers only the 
dynamics of sodium and potassium ion channels i.e. only these two channels are 
assumed to pose voltage dependent nonlinear conductance. Other ion channels are 
assumed to follow Ohm ’$ Law. The potassium channel is comprises of four activation 
gates whose opening probability is assumed to be n. The sodium channel has three 
activation gates with opening probability of m and a single deactivation gate with 
opening probability h. The circuit diagram of the model is shown in Figure 2.4. The 
dynamics of the model is given in Equation 2.7. 
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Figure 2.4: Circuit representation of Hodgkin-Huxley model 


C <1M± = G Na rn*h{E Na -V m (t)) + G K ri i (E K -V m (t)) + 
at 

Ch (El - V m ) + I in j O') (2-7) 

(1 - 8(t)) - h(VMt) ; a € n, m, h 
at 

a t (V m ) and p s {V m ) are voltage dependent rate constants. a s (V m ) determines the num- 
ber of transitions that occur from the closed to open states and (3 a (V m ) specifies the 
same from the open to closed states. These rate constants for the variables n , m, and 
h can be expressed as 




10 - v rn 


100 ( e ( 10 - V ' m )/10 _ 1 ) ’ 


Pn(V m ) = 0.125 e- Vm/S0 


LYui, ( ’On ) 


25 - V m 


10 (e( 25_v ' m) / 10 - 1) ’ 


0m(V m ) = 4 e'^ 18 


a h (V m ) = 0.07 e- vW2 °; Ph(V m ) = e(30 -v m )/ ~ io Tl 
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(a) Membrane potential (b) Time course of gating variables 

Figure 2.5: Response of Hodgkin-Huxley model 

The parameter values of HH model for a typical membrane are G Na = 120 mS/cm 2 , 
E Na = 115 mV, G k = 36 m,S/crn 2 , E K = -12 rnV,G L = 0.3 mS/cm . 2 , E L = 
0 mV [22]. Potential parameter values considered are taken with respect to the resting 
potential of membrane. Response of the model is plotted in Figure 2.5. Part (a) shows 
the variation of membrane potential with respect to time. The spike is initiated only 
when the current magnitude is more than some threshold value known as rheobase 
current. The time courses of m, h, n are shown in part (b). The sodium current is 
inward and depolarizes the membrane whereas the potassium current is outward and 
causes hyperpolarization. It has been observed that the response of n is similar to that 
of the membrane potential. The value of m jumps to its maximum when this potential 
crosses a threshold value. Incidentally, the value of h is also high at that time instant. 
This results in high magnitude of sodium current that leads to a high depolarization. 
Potassium channel is almost in deactivation state and could not hyperpolarizes the 
membrane. Therefore a spike is initiated. I he rise in membrane potential brings h, 
to a lower value and n to a higher value as seen from the time course of the gating 
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variables. The membrane potential hyperpolarizes owing to high outward potassium 
and low inward sodium current. The HH equations model both subthreshold and 
suprathreshold regime of neuron. This model shows the behavior of Type-II neuron as 
the firing rate is very high and spike latency is very low [5]. 

2.2.4 Morris-Lecar Model 

FitzHugh and Nagumo et. al. reduced the order of HH model and proposed a two 
dimensional model [12], [27] to make it computationally efficient. Although this model 
is studied for understanding the dynamics of membrane by several researchers, it has 
been argued that this model does not represent any practical neuronal membrane. In 
1981, Morris and Lecar suggested a two dimensional simplified model that describes 
the characteristics of barnacle muscle fiber [25]. This model incorporates two ionic 
currents, one outward going non-inactivating potassium current and the other inward 
going non-inactivating calcium current. Because calcium current responds much faster 
than potassium current, it is assumed that calcium current is always in equilibrium for 
the time scales considered. So, the dynamics of calcium currents can be ignored and 
the model is reduced to two dimensions. Following equations describe Morris-Lecar 
model: 

= -G Ca mM n ) ( V m (t ) - E Ca ) ~ G k -w(t) (V m (t) - E K ) - 
at 

G m (V m (t)-V ReM ) + I inj (t) (2.8) 

dw(t) _ WoojVm) ~ W(t) 
dl. r w (V m ) 

where w(t) is the recovery variable for potassium channel. The potassium activation 
variable follows the standard first order equation with steady-state activation 


WociVm) = 0.5 (l + tmh^tl 
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The time constant is 

MKn) = — ~ 777- 

cosh (Kn/60) 

It is assumed that the calcium current is always at equilibrium, with its activation 
curve given by 

'rnoc(Vm) = 0.5 ^1 + tank — 

There are many variants of the ML model [22] . However, the widely used one considers 
the following parameter set [35]: G Ca = 1.1, <2*- = 2.0, G m = 0.5 ,E Ca . - 100 ,E h - = 
—70, and V r ,.„ t = —50 (all conductances are in rnS/cm 2 and potentials in mV). The 




(a) Time course of variables (b) Phase portrait 

Figure 2.6: Response of Morris-Lecar model: Case 1 

response of the ML model is shown in Figure 2.6. Upper subplot of part (a) shows 
the firing characteristics. The shape of spike is very sluggish and interspike interval is 
high. This dynamics represents the properties of Type-I neuron. Lower subplot shows 
the time course of recovery variable w(t). A phase portrait between the membrane 
potential V m {t) and recovery variable w(t) is plotted in part (b). The phase portrait 
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Tim® (ms) v ( mV > 


(a) Membrane potential (b) Phase portrait 

Figure 2.7: Response of Morris-Lecar model: Case 2 

shows a limit cycle. To find out the effect of initial conditions on the ML neuron 
model, the analysis has een carried out for different initial conditions. The result is 
shown in Figure 2.7. As depicted from the Part(a) of the figure, there is a fail in 
generation of action potential (represented by the green line) at some initial condition 
(V£ = -21, w° = 0.09). Also, the parameter set allows only a single spike irrespective 
of initial conditions instead of repetitive spiking. The locus of the phase portrait is not 
a closed one unlike in the case of first parameter set. 

2.3 Synapse Modeling 

2.3.1 The Synapse 

The gap junction between the axon of one neuron and the dendrites of the next neuron is 
termed as synapse. The synapse allows the transmission of impulse from a neuron to the 
other. The synapse lies between two terminals called as presynaptic and postsynaptic. 





18 


terminals. Ihe axon ends at the presynaptic terminal and the dendrites begin from 
postsynaptic terminal. Synapse is classified as 

• Chemical Synapse 

• Electrical Synapse 

The space between the presynaptic and postsynaptic terminals of the chemical synapse 
is called the Synaptic Cleft.. A number of synaptic vesicles exist in the presynaptic 
terminal. These are spherical membrane bound organelles. Each of the synaptic vesicles 
contain one or more neurotransmitters. Neurotransmitters are secreted as a form of 
chemical signals from the presynaptic end. Neurotransmitters work as messengers 
that convey information from one neuron to the other. The neurotransmitters are 
classified into two categories small molecule neurotransmitters and neuropeptides. 
Small molecule neurotransmitters mediate rapid synaptic actions and neuropeptides 
tend to modulate slower, ongoing synaptic functions [46], 

Even though the chemical synapses outnumber its electrical counterpart, the elec- 
trical synapses are found in all nervous systems including human brain [31]. In case of 
electrical synapse, the membranes of the two communicating neurons come extremely 
close at the synapse and me actually linked together by an intercellular specialization 
called a gap junction. Gap junctions contain precisely aligned, paired channels in the 
membrane of the presynaptic and postsynaptic neurons so that each channel pair forms 
a pore. A variety of substances diffuse between the cytoplasm of the presynaptic and 
postsynaptic neurons. Electrical synapses thus work by allowing ionic current to flow 
passively through the gap junction pores from one neuron to another. The usual source 
of this current is the potential difference generated locally by the action potential. 

2.3.2 Dynamics of Synapse 

Synapse transmits the information through synaptic current. The value of current 
depends on the synaptic conductance g sy n(t ), synaptic potential E syn and the action 



potential V m (l.) arrives at presynaptie terminal. Synaptic current is governed by Equa- 
tion 2.9 


— 9syn{t) (V m (£) ~ E syn ) (2.9) 

The presynaptie potential is considered to follow the HH response as shown in Fig- 
ure 2.5. Chemical synapses a, re classified as NMDA Synapse and Non-NMDA Synapse 
[22]. NMDA is the abbreviated form of N-methyle-O-aspartic acAd. The synaptic con- 
ductance for the NMDA synapse is modeled using the Equation 2.10 

SsynO') ~ const • t e~ i,tpsak (2.10) 

tpeak is the time at which the synaptic conductance cj syn (t) reaches its maximum value 
The value of the constant is s ^~. For the values t p , uk = 0.5 ms and g pea k = 
1 nS, the variation in synaptic conductance and synaptic current is plotted in Figure 2.8 




(a) Synaptic conductance 


(b) Synaptic current 


Figure 2.8: Characteristics of NMDA synapse 


Major part of the synaptic current lies in the fourth quadrant. Hence, this synaptic 
current inhibits the next neuron. The value of synaptic conductance rises very fast. 
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The fall of conductance value is more flat when compared to the rising time. Synaptic 
current acts like a transient as it falls to zero within a very short period. 

The Non-NMDA synapse transmits the signal with the help of Mg 2+ ion. Synaptic 
conductance is given by, 

_ e -tjr 2 

9syn{t) = .< 7» 1 + ~ [ A / 9 2+] e - 7 V m (t) ( 2 - n ) 

The value of cj n varies from 0.2 nS to 0.4 nS at 35° Celsius. Physiological concen- 
tration of Mg 2+ ion is approximately 1 mM. The values of other parameters are 
T\ — 145.5 ms, T‘i = 4.1 ms, a - 0.33 /mM and 7 = 0.06 /mV. The synaptic current 
is obtained by directly multiplying the synaptic conductance to the membrane poten- 
tial as Esyn & 0 [22]. The time course for synaptic conductance and current are shown 
in Figure 2.9 




(a) Synaptic conductance (b) Synaptic current 

Figure 2.9: Characteristics of Non-NMDA synapse 

The Non-NMDA synaptic conductance is a function of action potential. It has a 
shape similar to a bell-shaped function. The current initially is inhibitory, then an 
excitatory pulse is generated, and finally it again goes to inhibitory state. However, 
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this current, acts as an excitatory current as the generated pulse is excitatory. The 
Figure 2.10 shows a current voltage relationship for Non-NMDA synapse. As the spike 



Figure 2.10: Current-voltage relationship for Non-NMDA synapse 


falls to the resting potential after it fires as shown in Figure 2.5, the current-voltage 
relationship shows two values of current for a particular membrane potential. The 
contour of the curve starts from point X. follows the path shown by arrows and ends 
at point Y. It intersects itself at Point A. Hence, the value of synaptic conductance is 
same for both rising and falling edge of spike at this membrane potential. At all othei 
cases, the synaptic conductance does not. take the same value for rising and falling edge 
for any potential value. 


2.4 Coupling in Neural Systems 

Determining the dynamical behavior of an ensemble of coupled neurons is an important 
problem in computational neuroscience. From the very beginning of the research m the 
field of computational neuroscience, people deal with single neuron and its behavior. 
Present trends of research include investigation of the behavior of neurons considered 
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in a network and their way to fire synchronously. It is assumed that the activities 
in the brain are synchronous and underlying interests for synchronization of nonlinear 
oscillators in physical and biological systems range from novel communication strategies 
to understand how large and small neural assemblies efficiently and sensitively achieve 
desired functional goal [30]. The coupling can be linear or nonlinear. The linearity of 
coupling on a neuron depends on the coupling current introduced by the other neurons 
coupled to it. Coupling can also be classified as strong and weak based on its strength. 
The impact of linear coupling on the neuronal behavior will be analyzed in the following 
sections. 


2.4.1 Linear Coupling 

The simplest coupling phenomenon is the linear coupling between two neurons. If 
Neuron A and Neuron D are coupled with each other, then a current will induced from 
Neuron A to Neuron B and vice versa. In case of linear coupling, at any time instant, 
this current depends directly on the difference between the potentials of the two neurons 
and the coupling conductance. The coupling conductance determines the strength of 
the coupling strength. In synaptically coupled neurons, the dynamics of one neuron 
affects the other only when the potential of that neuron crosses a threshold limit. 
Morris-Local’ model has been used to find out the coupling effect between neurons. 
The dynamics of the i th neuron in case of n mutually coupled neurons can be given by 
Equation 2.12 

- -5k< (K) («4<«) - EL) - Biyw (K(t) - ft) - 

o-. (km - ft) + t < 2 - 12 > 

j=l, 

dw i {t) _ wM)-^(yl) 

dt r 4(Kn) 
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The Heaviside Function , H (V$(f) - V Th ) is used to implement the threshold condition. 
The activation curve for the calcium and potassium conductance and the time constant 
used in Equation 2.12 are given by Equation 2.13. 


m’L =0.5 1 + tanli 


K+± 

15 


v i - 10 

uC. = 0.5(l + tanli-2— - 
14.5 


3 


cosh ha 


29 


(2.13) 


Table 2.1: Parameter sets for the neurons 


Neuron 

G Ca 

G k 

G M 

Eca 

Ek 

Em 

Neuron A 

1.2 

2.0 

0.5 

120 

-84 

-60 

Neuron D 

0.8 

2.3 

0.6 

115 

-77 

-65 

Neuron C 

1.0 

2.1 

0.4 

118 

-62 

-80 


It has been assumed that each of the coupled neuron follows Morris-Lecar model as 
seen from Equations 2.12 and 2.13. This dynamics is found to occur in specific part of 
brain mostly in the purkinje cell of cortical neuron. However, the parameter values for 
each neurons are kept different. Owing to nonidentical parameter values, each neuron 
poses different characteristics such as different firing frequency, rheobase current, and 
different refractory period. The differences in parameter values occurs due to different 
ion concentrations even if the neurons are closely located. For the present chapter, 
the dynamics of three neurons denoted by Neuron A, Neuron B and Neuron C lias 
been considered. The parameter values for the neurons are shown in Table 2.1. The 
membrane capaeitance is turned to be 1 pf/cm 2 for ail the three neurons. Therefore, 
the capacitance C is oninnted from the equations considered now onwards. 
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2.4.2 Dynamical Analysis 

In order to carry out the dynamical analysis, the equilibrium points are calculated for 
different stimulus magnitudes. It has been done by forcing the time derivatives of the 
state variables to zero. By setting w l = 0, the value of recovery variable at equilibrium 
point is obtained as 

.. ( V 1 - 10 \ 

w l = ■<, = 0.5 ( 1 + tanh - 

The other condition for the equilibrium is = 0. After substituting the value of 
w l * in this condition, Equation 2.14 is obtained which when solved gives the value of 
membrane potential at equilibrium point. 


-0.5 a Ct> 1 + tanh 


V ff + 1 

15 


(v£ - Ek) 

(v£ - K) - a,;, (q; - eL) + 


0.5G 


K 


V { * - 10 
1 + tanh ■ — — - ■■■ 
14.5 


E 9:?(vj 

; I, jrjLi 




■% 

inj 


0 (2.14) 


Equation 2.14 can not be solved explicitly for V£. A numerical approach has been 
taken to find out the values of for different values of I\ n y Once the locus of the 
equilibrium point with respect to current is obtained, the eigenvalues can be calculated 
to check the stability of the system. Eigenvalues are the roots of characteristic equation 


d,et\J-X I\ = 0 


where J is the Jacobian Matrix. Therefore, it is necessary to find out the Jacobian to 
evaluate the eigenvalues. The coupled ML system has the form: 
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The Jacobian is formed with the values of the derivatives of /' and cf with respect to 
its state variables calculated at its equilibrium point. Therefore, the Jacobian is given 
by, 



mi 

8V,}, 

Mi 

a(/2 

u 'm 

l£i 

8V* - 

Mi 

dw 1 

Mi 

dw 2 

M. 

dw' A 

jit 

a/ 2 
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Mi 
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Mi 
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Mi 
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This yields eight types 



are given by Equations 2.15 2.22. 


ML 

dw'J 


. These 


|ZL = -Q.5 (t Cll (l + tanh ) - 0.033 G' Ca ( 1 + 




(Vi - Bb„ ) - w‘ - Gi, - £ sf> II (vl - v n ) (2.15) 




d JJ = <?>H(Vl-V Tk ) 


(2.16) 




(2.17) 


df 


duJ 


T = 0 


(2.18) 


S = 0-012 ^ (°- 5 (* + tmh 1 - *"■ 


+0.012 cosh ( — — — ] I 1 + ( tanh 1 m 


10 


V 
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Vl - 10 


14.5 


(2.19) 
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dg‘ 

• =0 

OVA 

(2.20) 

dg i __ 
dw i 

-0.33 cosh 10 j 

(2.21) 


O 

II 

'*£>1 % 

(2.22) 


ihe Jacobian is formed by evaluating Equations 2.15 - 2.22 at the equilibrium points 
0 w ' w ‘ )• Eigenvalues from Jacobian matrix are calculated in order to determine the 
specific type of bifurcations in system. 

2.4.3 Two Neuron Coupling 

Consider the coupling between Neuron A and Neuron B . The independent response 
of the two neurons are drawn in Figure 2.11 when they are subjected under a stimulus 
of 20 nA and 35 nA respectively. 



Figure 2.11: Time course of Neuron A and Neuron B (no coupling) 

The red line shows the time course of Neuron A and the blue line represents the 
same of Neuron B. It is observed here that the interspike interval of Neuron A is much 
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(a) Time course of membrane potentials (b) Phase portrait 

Figure 2.12: Response of weakly coupled neurons 

larger when compared with that of Neuron B. The response for the coupled model with 
coupling conductance g 12 = g 2 } = 0.1 mS/cm 2 is plotted in Figure 2.12, It is found 
that Neuron A and Neuron B try to synchronize with each other as the difference 
in their firing frequencies reduces to 29.73 Hz from 33.96 Hz. The phase portrait 
between the membrane potential of each neuron is plotted in Part (b). To achieve 
synchronization, the firing frequency of the two neurons should be same and therefore 
the phase portrait between the membrane potentials must yield a straight line. An 
increase in coupling strength may drive the neurons into synchronism. The behavior 
of the neurons are studied for g 12 = gf = 0.8 mS/cm 2 and the results are shown in 
Figure 2.13. It is evident that the strongly coupled neurons are almost in synchronism. 
The two neurons fire together with a frequency of 121.5 H z. 

The variation in the behavior of the neurons has been observed by plotting bi- 
furcation diagram at different injected current. Two cases are considered, one with 
gi‘2 _ o.i mS/nm? for weak coupling and the other with g\ 2 = 0.8 mS/gnrv for strong 
coup lin g. The stimulus magnitudes are assumed to be identical for both the neurons, 
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(a) Time course of membrane potentials (b) Phase portrait 

Figure 2.13: Response of strongly coupled neurons 

i.e. I [ n;i = l™* — 1. The system has four system states and the state vector is given 

by, 

->f = fC V* m< .» ! f 

Owing to its four dimensionality, it has four eigenvalues for each equilibrium point. 
The equilibrium point shifts as the magnitude of stimulus is varied. This results in a 
change of the eigenvalues with respect to stimulus which helps to describe the system 
dynamics at the equilibrium point. Part(a) of Figure 2.14 plots the locus of eigen 
values with respect to stimulus. Two (Ai, Aa) of the four eigenvalues are pure real and 
the other two (A 3 , A4) are complex conjugates. The magnitude of real part of A x , and 
A 2 crosses the zero value for / = 12.95. The other two eigenvalues remain positive. 
This condition shows the existence of Hop} Bifurcation. The equilibrium points and 
eigenvalues for the corresponding stimulus are given in Fable 2.2 
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(a) Locus of eigenvalues (b) Bifurcation diagram 

Figure 2.14: Bifurcation analysis for two weakly coupled neurons 


Table 2.2: Values for weakly coupled neurons 


State 

Equilibrium Point 

Eigen Value 

V '1 

m 

7.74 

2.52 

V 2 

m 

34.46 

0.46 

w l 

0.43 

0 + j 2.78 

w 2 

0.97 

0 - j 2.78 
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The Jacobian at equilibrium point for the stimulus at bifurcation point is 

3.31 0.1 -183.48 0 

0.1 0.46 0 -256.37 

JgP-.O.l = 

0.01 0 - 0.33 0 

0 0.03 0 -0.46 

The bifurcation diagram portraits the change in system state when the bifurcation 
parameter is varied. The height at any particular point shows the change in magnitude 
between two time instants at that value of bifurcation parameter. Therefore, the system 
behavior changes with the sudden increase of height from a very low value. In the 
Part(b) of Figure 2.14, red dots show the bifurcation of Neuron A and the same for 
the Neuron B is shown by blue dots. There is a behavioral change, i.e. transition from 
stable to oscillatory mode at / = 12.95 nA. The result is a Hopf Bifurcation. However, 
the response of the Neuron B is not identical with that of Neuron A due to the weak 
coupling of neurons. Bifurcation Analysis is carried out for strongly coupled neurons 
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(a) Locus of eigenvalues (b) Bifurcation diagram 

Figure 2.15: Bifurcation analysis for two strongly coupled neurons 
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and the response is shown in Figure 2.15. The real parts of a complex conjugate 
eigenvalues cross the zero value for / = 6.31 nA. Hence, the Hopf Bifurcation occurs 
at this value of current. 1 he system moves from the stable state to oscillatory state. 
This is also evident from the bifurcation diagram. The red and blue dots show the 
responses for Neuron A and Neuron B respectively. The distribution of red and blue 
dots are equally dense as the system is synchronized because of strong coupling between 
them. Fable 2.3 shows the eigenvalues and the equilibrium points. The Jacobian for 
I = 6.31 nA is computed as 


JtjU-.M — 


2.35 0.8 -175.25 0 

0.8 0.65 0 -208.39 

0.01 0 -0.34 0 

0 0.01 0 -0.34 


Table 2.3: Values for strongly coupled neurons 


State 

Equilibrium Point 

Eigen Value 

V" 1 

y m 

3.63 

1.17 + j 0.52 

v 2 

Y m 

13.6 

1.17- j 0.52 

IV 1 

0.29 

0 + j 1.59 

w 2 

0.62 

0 - j 1.59 


2.4.4 Three Neuron Coupling 

In this section, the coupling behaviors of three neurons, Neuron A, Neuron B and Neu- 
ron C are analyzed in order to find out the effect of increased population on coupling. 
Stimulus applied to Neuron A and Neuron B are identical to earlier two neuron cou- 
pling model, i.e. 20 nA and 35 nA respectively. Neuron C is subjected under a current 
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of magnitude 10 n.A. I he parameters of Neuron C are chosen in a way such that it 
remains in the subthreshold regime for lower stimulus magnitude. The idea behind this 
is to find out the effect of coupling on neurons that do not fire under normal conditions. 
I he independent uncoupled response of the three neurons are shown in Figure 2.16. 
As mentioned before, the Neuron C could not cross the threshold and hence no spike 
is initiated in it (black line). 



Figure 2.16: Time course of Neuron A, Neuron B and Neuron C (no coupling) 


If the neurons are coupled together with a very low coupling conductance g l J = 
0.1 mS/cni 2 , the third neuron crosses the subthreshold regime and starts firing. How- 
ever, the whole system is highly unsynchronised as depicted from the randomness of 
the phase portrait shown in Figure 2.17. Time courses of the neurons also show dif- 
ferent firing frequencies. The Neuron B has highest firing rate followed by Neuron A. 
Although Neuron C exhibits firing, it has the lowest firing frequency. Therefore, the 
Neuron C needs more coupling current to fire synchronously with the other two. 

The coupling conductance for the three weakly coupled neuron model has the same 
magnitude of that of two weakly coupled model. The phase portrait between Neuron 
A and Ne.uron B is less random in the present condition when compared with the two 
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(a) Time course of membrane potentials (b) Phase portrait 

Figure 2.17: Response of three weakly coupled neurons 

neuron model. Similarly the firing frequency difference between Neuron B and Neuron 
A is slightly less (lesser by 1.57 Hz from two neuron case) in three neuron model. It 
can be concluded that the introduction of third neuron affects, though not in great 
extent, the synchronization between other two neurons. 

The next objective is to find out whether a higher coupling strength brings Neuron C 
into synchronism like the earlier two neuron coupling model. The coupling conductance 
has been increased to 0.8 mS/cm 2 . The neurons are now synchronized with each 
other. All the three neurons fire with a frequency of 109.82 Hz. The phase portrait 
between two neurons considered at a time is almost linear. Therefore the coupling also 
synchronizes the neurons that are uncapable to cross the subthreshold regime. The 
response is shown in Figure 2.18. 



Potential (mV) Membrane Potential (mV) 
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(a) Time course of membrane potentials (b) Phase portrait 

Figure 2.18: Response of three strongly coupled neurons 



(a) Weak Coupling 


(b) Strong Coupling 
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Like the two coupled noui on model, a bifurcation analysis is carried out for the 
present model. The stimulus has been considered as the bifurcation parameter. It is 
varied in steps to find out the changes in the behavior of the systems. A bifurcation 
diagram is drawn in Figure 2.19 with respect to stimulus magnitude. The red, blue 
and green dots show the membrane potential of Neuron A, Neuron B and Neuron C 
respectively. The coupled system enters in to oscillatory state from converging state 
at a current of 11.75 nA for weakly coupled and 24.8 nA for strongly coupled neurons. 
The diagram shows presence of Hop f Bifurcation. 


2.5 Summary and Conclusion 

In this chapter, the architecture and functions of a neuron and its different parts are 
described. The neurons are classified based on its topology, functioning and dynamical 
behavior. Characteristics of a passive membrane is studied. Response of Type-I and 
Type.- 1 1 neuron models are analyzed. The effect of initial condition is studied for Type-I 
neuron model. The Morris-Lecar neuron model shows a Type-I behavior whereas the 
Hodgkin-Huxley neuron model has a Type-II characteristics. The neuron models may 
exhibit a periodic spiking as well as a single spike. For any magnitude of stimulus, the 
neuron may and may not fire depending upon its initial conditions given by resting 
potential and the ion conductance at that potential. Synaptic response for both NMDA 
and Non-NMDA synapse are observed. For a presynaptic. potential that follows HH 
characteristics, the synaptic current is found to be inhibitory in case of NMDA synapse 
and excitatory in case of Non-NMDA synapse. Finally, the behavior when two and 
three neurons are coupled together has been examined. The magnitude of coupling 
strength affects the synchronization between the neurons. It is possible that a neuron 
may fire synchronously with other neurons even if it is uncapable to fire under lower 
stimulus magnitude. The coupling currents from other neurons help to fulfill the current 
deficiency required to fire the neuron. A bifurcation analysis has been carried out on 
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coupled neuron models to find out the effect of stimulus. The coupled neurons enter 
to oscillatory mode for the equal magnitude of current. However, in case of weakly 
coupled neurons, there exist a frequency difference between the firing of neurons unlike 
the strongly coupled neurons. The dynamical analysis show that the excitability of 
neurons are very much dependent upon the coupling strengths. The magnitude of 
coupling strength also shifts the bifurcation point. These dependencies have to be 
taken care while implementing the coupling phenomenon in real systems. 



Chapter 3 


Signal Propagation in Neuron 


The communication occurs in nervous system by the means of electrical signal. A 
neuron comprises of dendrites, soma and an axon. The dendrites receive electrical 
signals as a form of current from synapse. These current signals change the potential 
of dendritic membrane. Now the signals propagates down the dendrites as a form 
of voltage and summed up at soma. The aggregated potential is transferred to Axon 
hillock though a conductance. The current at Axon hillock initiates an Action Potential 
that travels along the axon which again converts into current signal at synapse that 
excites the dendrites of next neuron. Therefore, studying the phenomenon of signal 
propagation is very important in order to understand the behavior of nervous system. 


3.1 Cable Theory 

3.1.1 Cable Equation 

The dendrites and axon have a shape similar to a tapered cylinder. They are referred 
as cable owing to their cylinder-like structure. Since the diameter of this cylinder is 
very small, there is physically no variation in potential along the cross-section of it. 
However, the potential gets delayed in both active and passive membrane and decayed 
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in case of passive membrane as it propagates. The delay depends on its propagation 
speed. The decay in passive membrane occurs due to the longitudinal resistance of 
the cable. The passive membrane does not have any ion channels embedded in it that 
prevents the attenuation. The working of axon and dendrites are explained with the 
help of cable theory. The cable equation is given by, 


,W. i d ( 2 dv m (x,t) 


dt 


2(177, dx 


dx 


i"in (*E 1 1 ) "t" I inj (*r 1 1 ) 


(3.1) 


where r t . is the intracellular resistivity and a is the radius of cable. The membrane 
current may follow either the Hodgkin-Huxley dynamics or simple passive membrane 
model depending upon the part of neurons considered. The Equation 3.1 is reduced to 
Equation 3.2 when the cable has uniform cross-section. 


„c)V„,{x,t) a d 2 V m (x,t) 


' 'hn(* i 0 "b I'inj (•£, t) 


(3.2) 


()l 2 r L dx 2 

In case of passive membrane, the membrane current is given by i m (x , t ) = V m (x, t ) / R , n , 
where R, m is meinbrane resistance. This reduces the cable equation to 


CIl, 


dV m (xk) aR m d 2 V m (x,t) 


V m (x,t) Ij n j (.X' , t) R 


(3.3) 


dt 2 r L dx 2 

The potential is a function of time and space. The passive membrane does not 
have any ion pumps embedded in it that prevents the decay in potential. Therefore 
the passive membrane potential is attenuated both with respect to time and space. 
This potential lias a time constant r = CR. m and space constant £ = yj^. A typical 
membrane has a time constant r = 1 — 2 ms [22] and space constant £ ~ 1 mm [9]. The 
signal is offered a resistance when it travels longitudinally. This longitudinal resistance 
is given by /f £ = ^ . The analytical solution for the passive membrane model 

may lie found out for some current waveforms. 
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1 he < oust ant, ( un ent injection of magnitude I e yields a solution given by, 

V, u (.r) = 

1 lu i active meml uanos like axon that, have ion pumps embedded to facilitate the 
propagation can not be solved analytically as the dynamics of ion channels introduce 
high nonlinearity in it. Hence a numerical approach (multicompartment model) is used 
to model axon or dendrites. 

3.1.2 Multicompartment Model 

This is a numerical method to compute the action potential. Instead of considering the 
axon or dendrite as a. whole, they are split in to many segments. Each segment has a 
very small length and therefore assumed to be equipotential along its length. In other 
words, t he potential does not change along the length of each compartment considered 
separately. Two successive compartments are coupled with a conductance and hence 
t he potential differs due to drop in the conductance. The ion pumps modeled by HH 
equations inject, current, in order to prevent the decay in action potential so that it 
can propagates down to the presynaptic terminal. The current equation for the g t/! 
compartment is given by Equation 3.4 

C d }jp. = -4(f) + <r*-' (Kr'W - TO) +9 h * 1 - V£(t)) + 

W (3.4) 

The coupling conductance between two successive compartments g and g' is 

^ 

J r L L fl {U l a^ + />V 2 ) 

where a>‘ is radius and // is the length of the g th compartment. Owing to very 
small length of each compartment, the variation in potential across the length of com- 
partment is negligible. Thus the spatial derivative of the potential vanishes from Equa- 
tion 3.4. The partial differential equation resolves into N coupled nonlinear differential 
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quations. In order to find out the solution of Equation 3.4, the potential of the {^.+l) th 
ompartmeut at the same time instant must be known. Therefore, an implicit method 
ueh as Backward Euler or Crank Nicholson Method may be used. The first method 
< more stable whereas the second is more accurate. The results are shown from the 
nalysis omicd out with C rank Nicholson Method. The Backward Euler method and 
7 rank Nicholson method are described in Appendix-A. 


5.2 Dendritic Propagation 

"he dendrites do not contain any ion pumps. They tire modeled using a capacitor 
n parallel with a resistance. Owing to its passive nature, the signal attenuates as 
t propagates toward soma. The dendrites are very small in length and larger in 
ross-seetional area. Hence the attenuation does not result in complete die out of 
ignal. The pot ential reached at soma from different dendrite terminals are aggregated 
ogether that, stimulate the Axon hillock which in turn initiate the action potential, 
.'he dynamics of the dendrite is governed by Equation 3.5 

a &V n (x,t) V m {x,t)-V R e at , r 

( ~ 2r h ~ n? Wn + inj{ ’ ] [ ] 

The break up of dendrites into small compartments results in to the change of 
he Equation 3.5 to Equation 3.G for the /.i th compartment. The synaptic current is 
mlsative in nature as observed from Figure 2.8 and 2.9. The response of a dendrite 
objected under a pulse current is shown in Figure 3.1. The major part of the den- 
lrit.es lit's in synaptic area. This helps stimulus not only to be injected at the dendritic 
erminal, but also to spread along its stem. This in turn reduces the decay of signal. In 
he simulation, different fractions of stimulus are injected from first to third compart- 
tient. The potentials show sharp peaks for these compartments as seen from Figure 3.1 
iiid for the later compartments, the peaks are smoother. If the dendrites would have 
>eeii stimulated only at the first compartment, the decay results in complete die out 
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of potential. 
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(a) At different locations with respect to time (b) Three dimensional plot 


Figure 3.1: Signal propagation in dendrite 


3.3 Axonal Propagation 

Tht' signal in membrane are of two types Graded Potential and Action Potential. 
Graded potentials are short lived depolarizations or hyperpolarizations of an area of 
membrane. These changes cause local flows of current that, decrease with distance. 
Axon is very long in length and hence the graded signal die out completely if it travels 
down the axon. To transmit the information from a neuron to other, the Axon hillock 
generates action potential. The generation and transmission of action potential is 
facilitated by the ion pumps embedded in it. 
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3.3.1 Propagation of Action Potential 

In order to find out the behavior of action potential when it travels down the axon, 
a. multieompartmental analysis is carried out. The axon that has been considered for 
analysis has a length of 12 cm and a tapering cross-section with the diameter varies 
from 8 fim near soma to 2 fin: at its end. It has been split into 25 compartments. The 
action potential 1 ( 0 at. // compartment follows Equation 3.4. The result is shown 
in Figure 3.2 




0 0 


(ft) At different locations with respect to time (b) Three dimensional plot 
Figure 3.2: Action potential propagation in axon 

As observed from the figure, the spikes at the distant locations occur at a later time 
instant. Thus the action potential is delayed and occur at a later time. The shape and 
magnitude of the potential is identical except small perturbations like slight increase 
or decrease in magnitudes, more or less sharpness etc. The uneven ion concentrations 
in different ion channels are responsible for these perturbations. 
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3.3.2 Noise in Axon 

Und.T normal eondit ions, a minimum value of current is required to generate an action 
potential. For the topology of axon considered, an action potential is fired for /P . = 

inj 

11 n.\ and it travels down the axon without being attenuated. If the current injection 
at Axon hillock reduces to 10 nA, the action potential fails to generate as shown in 
Figure 3.3. However, there may still be a chance that the action potential is generated if 
it suffers by noise. Noise may be a result of either of the two or both reasons mentioned 
below. 

• Uneven opening and closing of gates of ion channels embedded on axon 

• Epileptic Interaction [18] by a nearby axon. 

1 lu* noise is assumed to follow a Gausnian distribution with zero mean and unity 
variance />. a white gaussian noise. The stochastic multicompartment model of axon 
is given by Equat ion 3.7 

0 +/'"" 1 (Kn*) “ VM)) +<f" l+l - V£{tj) + 

+ (3-7) 

I? nj (t)= 0 i/PTO 

The probability of occurring action potential has been found out for a range of 
noise intensity with power output varying from 1 dBW to 5 d,BW. The results are 
shown in Table 3.1. The part(b) of the Figure 3.3 shows the generation of action 
potential at / - 10 nA with 4 dBW noise power. It is clear- from the table that noise 
does not affect t he generat ion of action potential to a great extent. The 5 dBW white 
gaussian noise initiates an action potential if the current at axon hillock is just lesser 
than the threshold current required for generation. However, it fails to do the same if 
the current drops to 9 nA. Therefore the system is very less prone to noise and the 
system is robust. There would not Ire a false firing when the signal is lower than the 


threshold condition. 



44 


3. 1 : l-Jfcct of noise power on generation of action potential 


Samples 

Noist 1 Po war 

Spike Initiation 

%Spiking 


(ilBW) 

I,nj - 10 nA 

— 

II 

I in j — 10 Tl A. 

hnj 

= 9 nA 


1 

18 

0 

36 


0 


2 

15 

0 

30 


0 

5( 1 

[ 

3 

24 

0 

48 


0 


1 

2!) 

4 

58 


8 


o 

A .. • ■ 

30 

3 

78 


6 




(a) No noise (b) With 4 dBW noise power 


Figure 3.3: Action potential at / = 10 nA 





3.4 Fourier Analysis 
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The inullh ompmtmcut model is wry complex. It requires implicit method for analysis 
and the stability o! the solution depends on many factors viz. integration time step, 
number of < ompait moots and type of method. The shape of the action potential does 
not change with time. The propagation introduces only delay. It has been considered 
in spiking neuron model that the shape of action potential does not contain any infor- 
mation. Hie drawbacks of multicompartment model along with this phenomenon is 
tin* inspiration of Fount r Analysis. The Fourier analysis can be performed only in pe- 
riodk waveforms. The spike-ditfuse-spike model, which is observed numerically in more 
detailed biophysical models, has periodic traveling wave solutions [8]. The Figure 3.2 
also shows the periodic nature of action potential. Therefore the action potential can 
bo modeled using Fourier series. The point neuron HH model has been decomposed 
into a series combination of sinusoidal terms. These terms are modified with proper 
transformation in order to achieve the waveform similar to that of an action potential. 

3.4.1 Fourier Representation of Hodgkin-Huxley Neuron Model 

The dynamics of HH model is governed by Equation 2.7. Simulations have been carried 
out for different values of injected current ranging from 10 nA to 40 nA. It is found 
that the frequency of the spiking increases with the increment in the magnitude of the 
injected current. As the frequency plays an important role in the decomposition of any 
periodic waveform into Fourier components, an estimation of time period for different 
current injection is necessary. It has been found that the repetitive spiking occurs at a 
current value larger than 9 nA. A parabolic curve approximates the interspike interval 
(ISI) for different current values. The curve suggested here is given by, 


IS! = ( 13 - v /l.32(/„ 1J (f) - 9) )H(Iinj(t) - 9) (3.8) 

Heaviside function is used to take care that the Fourier series method is applied 
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Figure 3.4: Interspike interval as a function of injected current 


only in case of repetitive spiking. Figure 3.4 demonstrates that the curve given by 
Equation 3.8 can be used as an approximator of the time period. The blue line in the 
figure represents the actual curve of the time period for different current injection in 
HH model and the trajectory of the parabolic curve (approximated curve) is shown by 
the red line. 


The Fourier series of a function /(f) with time period 21 is given by Equation 


fit) 


CL 0 


Jforf 


E /Wl 1/ V N , 

a k cos— + 2^ b k stn 


font 


k= 1 


fc=i 


3.9 

(3.9) 


where, 

I r 21 

«o = y / f(t)dt\ 
l Jo 

(Ik =7 f f(t)cos^~dx; (3.10) 

l J o I 

1 /-2i ,, kid , 

bk== l Jo fi^ sm ~r dx 


Equation 3.10 considers the case when the function value is known in the domain 
0 to 21. The HH equations are coupled differential equations and hence the function 
/(f) is not known explicitly. However, the value of the membrane potential at different 
intervals can be obtained by solving the HH model numerically. Hence, to find out 
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values of oq, and b the integrals have to be subdivided into many intervals for 
which the value of membrane potential is known. Each of the integrals can be integrated 
numerically and added together to get the value of the coefficients. The Simpson's 1/3 
rule is used for the analysis as it is very simple and also gives fairly good results. 



Tim# (rna) Time (ms) 


(a) Individual Fourier components (b) After aggregation 

Figure 3.5: Membrane potential with Fourier analysis 

Figure 3.5 shows the response after the Fourier analysis for a current injection 
of 10 nA. With the help of numerically calculated Fourier components, the coupled 
nonlinear differential equations are represented algebraically. 

3.4.2 Constructing Action Potential 

Equation 3.1 is a partial differential equation in which the dependent variable action 
potential is a function of both time and space. The equation represents a wave. 
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A wave that travels at a. speed c in the x direction may be represented as 


V m {x, /.) = V m (ct - x) (3.11) 

The Fourier coefficients are obtained for the point neuron model. These coefficients 
can be used in a tactical way to compose a wave. The proposed model to generate the 
action potential is given by Equation 3.12 


,. 00 OG 

V m (x, t.) = -“ + 53 a kCos ( k ( ct - x)) + 53 hsin (k ( ct - x)) (3.12) 

Z k = 1 k = 1 

Equation 3.12 introduces the delay depending upon the speed of the propagation 
as the action potential moves along the axon. There is no attenuation and change of 
shape. It is solved numerically to obtain the final response of the propagating signal. 
The response is shown in Figure 3.6 




0 0 


(a) With respect to time at different locations (b) Three dimensional plot 


Figure 3.6: Action potential with Fourier analysis 
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The response is similar to that of multicompartmental model shown in Figure 3.2. 
However, there are some drawbacks in this Fourier series method. These are — 

• The shape of the spike slightly detoriates. It becomes more flat unlike the sharp 
spike of action potential obtained from the multicompartment model. 

• The response shows an existence of spike at axon end before the occurrence of 
spike near axon hillock. This has to be eliminated in order of proper functioning 
of the model. 

'1'he advantage of this model is that it is very simple and can be used robustly if 
the initial spike at axon end is suppressed. 

3.5 Analysis on Complete Neuron Structure 

The researches in Computational Neuroscience are carried out on point neuron model, 
dendrites, axon and synapse separately. The focus on the complete structure of the 
neuron has till now got very less significance. The dynamics of complete neuron is very 
much important for population analysis. 

3.5.1 Structure of model 

The neuron is considered to have three dendrites - Dendrite A, Dendrite B, and Den- 
drite C. Dendrite A and Dendrite B are assumed to receive the signal from Synapse 
A whereas Dendrite C receives signal from Synapse B. Synapse A is a Non NMDA 
Synapse and Synapse B is an NMDA Synapse. Furthermore, it has also been assumed 
that the dendrites are closely located and therefore no significant resistance exists be- 
tween two dendrites. The signals from synapse propagate through dendrites and are 
summed up at soma. This generates an action potential at Axon hillock. The axon is 
considered to be unbranched. 
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Two cases have been considered when the axon terminates at synapse. Analysis on 
the model has been carried out for each part separately. The cases are - 

• Termination at Synapse C with characteristics similar to that of Synapse A. 

• Termination at Synapse D with characteristics similar to that of Synapse. B. 

3.5.2 Dendrites 

The dendrites are modeled as a passive membrane. As discussed earlier, the dendrites 
take inputs from two synapses. The waveform of synaptic current is shown in Figure 3.7 




(a) Synapse A (b) Synapse B 

Figure 3.7: Synaptic inputs to Dendrites 

A nerve cell generally has a huge number of dendrites. It is not necessary that a 
neuron receives signal only from one synapse. In order to maintain the generality, two 
synapses have been considered that pose different characteristics. The topology of the 
dendrites are given in Table 3.2. 

The responses of Dendrite A and Dendrite B are shown in Figure 3.8. The response 
of Dendrite A is very much different from Figure 3.1. Two factors are involved for 
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Table 3.2: Topology of different dendrites 


Part 

Input 

Length 

Diameter near synapse 

Diameter near soma 

Dendrite A 

Synapse A 

1.0 mm 

12 pm 

16 pm 

Dendrite B 

Synapse A 

1.2 mm 

8 pm 

13 pm 

Dendrite C 

Synapse B 

0.8 mm 

10 pm 

12 pm 




(a) Dendrite A (b) Dendrite B 

Figure 3.8: Potential propagation in Dendrite A and Dendrite B 


wife m A,X£jL(J£'2m ~- —«* 
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the change in response of the dendrite the synaptic current and the time constant. 
Synaptic current is not exactly pulse shaped as considered in earlier case and also 
repetitive. Moreover, the time constant of the passive membrane is much larger than 
that of the synaptic current. This does not allow the dendritic potential to fall much. 
Also, there is no significant attenuation in the potential as it propagates. The synaptic 
current is negative for a larger time instant. It is distributed along the length of the 
dendrite instead of the injection only at its terminal. The result is almost similar 
potential waveform irrespective of space. 

Response of Dendrite B shows a clear rising and falling edge. The synaptic current 
distribution along its length is somewhat different from that of Dendrite A. The time 
constant is also lesser in case of Dendrite. B. The attenuation or growth in potential 
with respect to space is more significant. The reason for growth in potential is same 
with earlier case. Waveshape of potential propagation in Dendrite A and Dendrite B 
are similar with the difference only in rising and falling time. 



Figure 3.9: Potential propagation in Dendrite C 
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Dendrite. C receives the signal from Synapse B which is NMDA synapse. The 
response is shown in Figure 3.9. The potential near soma falls below the resting 
potential when the time is 8.14 ms. It remains negative with respect to the resting 
value up to t = 13.96 ms. This is due to the high inhibitory synaptic current between 
2.6 ms to 13.8 ms. 

3.5.3 Soma 

The signals from dendrites reach to soma as a form of potential. The absolute magni- 
tudes of these potentials are negative with respect to zero reference voltage. However, 
the effective magnitudes i.e. the magnitudes with respect to the resting potential may 
be positive or negative. The effective potential near soma for i th dendrite is given by, 

‘yjjcfj active __ -y jjabsolutt __ ^ 



Figure 3.10: The model of soma 

A schematic diagram showing the soma model is given in Figure 3.10. The voltage 
sources VD\, V D 2 and VD 3 are the effective potentials of Dendrite A, Dendrite B 
and Dendrite C respectively near soma. RD\, RD 2 and RD$ are the corresponding 
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Figure 3.11: Current injection at axon hillock 


resistances through which the dendrites are connected to soma. The summed potential 
Vs at soma appears across the resistance R. s . The soma acts as a buffer circuit that 
not only sums up the dendritic potentials and passes through the axon but also isolates 
the axon from dendrites. The current at Axon hillock generates the action potential. 
T his current shown by KVs is generated due to the multiplication of summed poten- 
tial Vg with the conductance that connects soma with Axon hillock. The current at 
Axon hillock is computed by carrying out the nodal analysis. This current is given by 
Equation 3.13 


I AH = K 


YJ2i + ¥J2z + 

RDr T RD2 _ 


VPs 

RDs 


1 

RD i 


+ ik + ^ + 1 


(3.13) 


RD 3 ' Rs 

where K is a constant that depends on conducatnce between soma and Axon hillock. 
The time course of this current is plotted in Figure 3.11. The magnitude of the current 
increases with the time from 0 ms to 3.7 ms. Then its starts decreasing up to 11.8 ms 
as the effective potential of Dendrite C is negative for this time span. The effect of 
this current on the generation of action potential is discussed in the next section. 
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3.5.4 Axon 

The axon considered for the analysis has a length of 12 cm with it diameter tapers from 
8 fim near soma to 2 /.mi at its end near synapse. The axon has no branches spreading 
out of it. Multicompartmental method is used for modeling the axon. Response is 
shown in Figure 3.12. It has been found earlier that there was no initiation of spike 
in ease of continuous current injection until the current magnitude exceeds 10 nA. In 
this case, the maximum magnitude of the current reaches 3.7 nA within the analysis 
time span. The topology of axon is also same in this case as in earlier case. However, 
the action potential is generated and propagated for the complete neuron model. The 
delay is also reduced. For the first spike, there is almost no delay when it propagates 
from 4 cm to 12 cm. The only change is in the shape of injected current. Therefore 
it can be concluded that the action potential may be generated for a lesser intensity 
of current if its magnitude is rising with respect to time. The action potential that 
arrives at the synapse is allowed to participate in the synaptic activities in order to 
compare it with earlier synaptic response. 



Figure 3.12: Response of axon in complete neuron model 
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The responses for the synapses are shown in Figure 3.13. As mentioned earlier, two 
eases have been considered. Part (a) of the figure shows the response when the axon 
is terminated at Non-NMDA Synapse and Part(b) shows the same when it terminates 
at NMDA Synapse. The characteristics are found to be almost similar except the 
frequency of occurrence of the synaptic current. 




(a) Synapse C (b) Synapse D 

Figure 3.13: Synaptic current 

A three dimensional plot for the potential distribution with respect to time and 
space for the three dendrites and axon are shown in Figure 3.14. 
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(a) Dendrite C 


(a) Dendrite A 


0.6 

0.4 

Distance (mrr$- 2 
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(b) Dendrite B 


(b) Axon 


Figure 3.14. Three dimensional plot for potential distribution in dendrites and axon 
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3.6 Summary and Conclusion 

In this chapter, the behavior of neuron is analyzed when the signal propagates in 
it. The cable theory has been used to model the dendrites and axon. The graded 
potential propagation in a passive dendritic tree under constant current stimulation 
is analyzed. It has been considered that the stimulus is not only applied to its end, 
rather it is distributed at its extensions near the synapse. The effect of stimulus with 
constant magnitude is also studied to model the active membrane of the axon. The 
study includes the effect of noise. A new method based on Fourier Series is proposed 
to model the axon. Finally the behavior of the complete realistic neuron model has 
been investigated. 

It lias been found that the behavior of propagating signal is highly dependent upon 
the nature of membrane. There is no initiation of spike in case of passive membrane. 
The passive membrane not only introduce the delay in signal, but also attenuation. On 
the other hand, there is occurrence of spike when signal propagates in active membrane. 
Although a delay occurs depending upon its speed, there is no attenuation observed. 
For the considered topology of axon, a steady current of magnitude 11 n,A initiate 
action potential. A lower value of current magnitude may also result in spiking if noise 
is induced. However, the model of axon is very robust or in other words immune to noise 
as a noise with 5 dBW power is also almost incapable to generate action potential when 
the current is 9 nA. The initiation of action potential also depends on the shape of the 
stimulus waveform. The Fourier Series method may be used to avoid the complexity 
of multicompartmental model for approximated work. The complete neuron modeling 
provides the details of information processing in neuron. The waveform of Synapse A is 
identical with that of Synapse C . Similarly, the waveforms of Synapse B and Synapse 
D are identical. The difference is only in the magnitudes and in time of occurring 
of current pulses. The Synapse A and Synapse B considers the dynamics of a point 
neuron HH model as a presynaptic potential and hence is the difference. The identical 
response of the synapses cross checks the validation of the model. 



Chapter 4 


Conclusion 


4.1 Work Done 

The neurons are coupled together to work synchronously for performing a specific task. 
In this work, the effect of coupling strengths on the network of neurons is studied and 
the following conclusions are drawn - 

• The weakly coupled neurons are not in synchronism as the coupling current 
that being proportional to potential difference between the neurons, controls the 
synchronous behavior of the neurons, is less. 

• The coupled neurons are in complete synchronism for strongly coupled system 
and there is no phase difference between successive firing of the neurons. 

• The variation in coupling strength results in different firing frequency of both 
the neurons. The firing frequency is more in case of strongly coupled system 
compared to weakly coupled system. 

• Synchronization phenomenon and variation in excitability is depicted by plotting 
bifurcation diagram for the pair of coupled neurons. 


59 



A similar analysis has also been carried out for the interconnections of three neurons 


and the following observations are made according to simulation results. 

• In ease of weakly coupled system, all the three neurons show spiking behavior 
whereas a neuron was uncapable to cross the subthreshold regime when not cou- 
pled with the other two neurons. The coupling mechanism provides additional 
stimulus which brings the neuron to oscillatory from stable state. 

• Similar to two coupled neuron case, the three strongly coupled neurons fire syn- 
chronously. 

• Bifurcation analysis is carried out on coupled system and it is found that for 
weakly coupled system with variation in bifurcation parameter, the neurons enter 
into oscillatory state earlier as compared to strongly coupled system. 

In Chapter 3, signal propagation in neuron structure is studied and different ap- 
proach for investigating the signal propagation is presented. The following conclusions 
are drawn from this analysis. 

• The addition of white gaussian noise reflects the biophysical activities such as 
ephaptic interaction and uneven opening and closing of ion channels and this 
facilitates in signal transmission by slightly reducing the threshold. Thus noise 
helps in firing, but not in great extent. 

• A Fourier series based methodology for signal propagation in neuron is proposed. 
This reduces the computational complexity. 

• In last section of this chapter, an overall signal propagation in a neuron structure 
by considering the dynamics of all the components are presented. It is found that 
the current at Axon hillock is not steady and variation in current at Axon hillock 
results in variation in spike speed. 
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4.2 Future Scope 

• The behavioral analysis on the coupled neuron model can be extended to find 
out the existence and effect of Chaos for a large population of neurons. 

• Study can be carried out on the synchronization of nonlinearly coupled neurons. 

• In the complete neuron model, the dendrites me assumed to be located very 
closely and therefore the resistance between two dendrites are ignored. The 
study of this model with the inclusion of these resistances may capture some new 
characteristics. 

• The effect of dendritic spines on the response of complete neuron can also be 
investigated. 

• As the dynamics of a complete neuron model is revealed in this thesis work, the 
analysis can be carried out on a population of neurons. 
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Appendix A: Implicit Integration Method 


The numerical solut ion of multicompartment model requires the value of potential at a 
later time stage. This requires implicit integration method. Two types of implicit inte- 
gration method are widely uned- Crank- Nicholson method and Backward Euler method. 
The membrane potential for the ji ,h compartment, is given by 

= -4(t) (V Z~\t) - V£(t)) + g >^ +1 iyT\t) - V£(t)) + 

Hnjit) 


This equation can be written as 


dt 


A ll V>r l + B ,L V!n + CNV£ +1 + D> 1 


where 


= g ^~ l / C ; = ($*'*" 1 + £>'' + g Ml+l ) fC 


Cf = g^ +1 /C; D>‘ = (</.?£* + lQ /C 

This generates a set of coupled differential equations. For z € [— 1, 1] 

V£{t + zAt)**V£(t) + z&V> 

The finite difference method yields, 

A VZ(t.) = + zAt) + B*V£[t + zAt) + CNV£ +l (t + zAt ) + If) At 

Approximating, 

A V£ = aPV£- 1 + b fl V£ + c?V£ +l + # 

where 

at 1 = A^zAt\ If = B^zAt; <f = C^zAt 
d fl = (D> L + AfV^t) + B»V£(t) + CSV£ + \t)) At 
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For N compartment, model, a 1 = 0 and c N = 0. The potential change for the (p - l) th 
compartment, is given by, 


, d 1 " 1 AV' 1 + & L ~ l 

a r//*-l _ c 2Z * m Z u new 
m I _ 

1 OtUiUJ 


where, 


i,e+i = ;/'+i + .—.Iff.,, ■ 
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r j/‘+ l = f /e+l . fl/ + ^new , ^1 ^1 

“new w ' “ i u* 5 u ncw u 

1 ~ Onew 


At the end of the cable, 


r1 N 

a yN _ u 'new 
^ Vm 1 _ AiV 


The algorithm for solving the multicompartment model is: compute bfi ew and d% ew 
down the length of the cable, then solve for AV^ and iterate back up the cable. This 
process takes 2 N steps. 

For Crank-Nicholson method , z — 0.5 and for Backward Euler method, z — 1.0. 
The Backward Euler method is more stable whereas the Crank-Nicholson method gives 
more accurate result. 



